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ABSTRACT 

We use cosmological hydrodynamic simulations to investigate formation of galactic 
' bulges within the framework of hierarchical clustering in a representative CDM 

o ' 

5-H 



cosmological model. We show that largest objects forming at cosmological redshifts 
z ~ 4 resemble observed bulges of spiral galaxies or moderate size ellipticals in their 



Q-e general properties like sizes, shapes, and density profiles. This is consistent with 

observational data indicating the existence of "old" bulges and ellipticals at more 
moderate redshifts. These bulges are gas dominated at redshift z = 3, with high rates 
! of star formation and would appear to be good candidates for small blue galaxies seen 

in the Hubble Deep Field. 

m 
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Subject headings: cosmology: theory - cosmology: large-scale structure of universe 
' galaxies: formation - galaxies: intergalactic medium 



Introduction 



Pl,: 

6 

' The fact that bulges of ordinary galaxies (and indistinguishable ellipticals of the same 

luminosity) are very dense is often used as an argument against the currently favored CDM-type 
cosmological models. Really, since the average density of the universe decreases with time, and the 
^ . average density of a bound object is directly proportional to the density of the universe at the time 

when the object is formed, dense galactic bulges should have formed at very high redshift. Thus, 
it is reasonable to ask whether currently fashionable cosmological models normally have enough 
small scale power to account for the formation of massive (10 9 - 10 10 solar masses in baryons) 
bulges at z £ 10 (Peebles 1997). 

This argument can be illustrated by the following simple estimate: the characteristic number 
density of baryons in the Galactic bulge within the sphere with the radius of 3 kpc is 

riGB ~ 2 cm -3 (1) 
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(Kent, Dame, & Fazio 1991; Freudenreich 1998). In comparison, the average density of the 
cosmological virialized object formed at redshift z is only 

ra T H ~ 10 _4 (1 + zfcm.- z (2) 

for Q b h 2 = 0.02, and we assume that the average density of the virialized object is about 200 times 
the average density of the universe at the moment of formation (Gunn & Gott 1972), according to 
the standard dissipationless collapse theory. Comparing equations @) and (H) we can deduce that 
the bulge of our Galaxy formed at zqb ~ 30. But, none of the currently acceptable CDM-type 
models can form 10 10 solar mass baryonic objects at z = 30 in numbers even closely comparable 
to the observed number density of galactic bulges. 

Does this argument imply that the CDM-type models are ruled out? We will try to show in 
this paper that the answer to this question is no. What this simple argument misses is the ability 
of baryons to cool and collapse to the densities exceeding that of the dark matter. We present 
a series of cosmological hydrodynamic simulations which include adequate physical modelling to 
properly account for formation of bulges of galaxies, and we show that typical bulges form in a 
realistic CDM-type model at z ~ 5 rather than at z ~ 30. 

But before describing our detailed results, it is perhaps worth noting the conceptual flaw 
in the simple argument we first presented. The global ratio of baryons to dark matter in our 
simulation, and in typical current estimates is about 1:8. But in our own galactic bulge the 
baryonic - stellar - component exceeds the dark matter component and may exceed it by as much 
as a factor of a few. Thus, (taking the cube root of density enhancement) the baryonic stellar 
component has contracted relative to the dark matter component by a factor of two to four and, 
correspondingly, the bulge should have formed in the plausible redshift range of z ~ 5 — 10 rather 
than at z = 30. Furthermore, as we shall see, the observed (and computed) profiles are very 
different from those envisioned in the top-hat collapse picture. 

An aside on a related but quite different problem may be useful here. There is currently a 
good deal of discussion on the subject of whether or not standard CDM models make bulges which 
are too dense (e.g. Spergel & Steinhardt 1999, Burkert & Silk 1999, Moore et al. 1999, Kravtsov et 
al. 1998). However, these papers address far lower mass systems and ones which are dark matter 
dominated, having circular velocities of less than lOOkm/s. Here we address normal bulges with 
(equivalent) circular velocities of about 200 km/ s or more, for which the advertised problem has 
been the opposite: why are they so dense? 

Some work on the subject was done by Steinmetz & Mueller (1994, 1995). However, they 
modeled the formation of disk galaxies as a collapse of an isolated gas cloud, whereas in this 
work we consider the formation of galactic bulges within the framework of hierarchical clustering, 
based on the realistic cosmological simulations that include such effects as the cosmological infall, 
merging of proto-galactic clumps, and expansion of the universe, all of which are missing in the 
Steinmetz & Mueller work. Other numerous numerical investigations of galaxy formation do not 
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specifically concentrate on the bulges, or lack mass and/or spatial resolution to say anything about 
bulge formation. 

2. Simulations 

We use the SLH cosmological hydrodynamic code (Gnedin 1995, 1996; Gnedin & Bertschinger 
1996). Physical modelling included in the code is fully described in Gnedin & Ostriker (1997). We 
choose a CDM+A cosmological model with the following cosmological parameters: 

Q = 0.37, Q L = 0.63, h = 0.70, n b = 0.049, 

which is close to the "concordance" model of Ostriker & Steinhardt (1995) and to the models 
consistent with recent SNIa results (Riess et al. 1998, Perlmutter et al. 1999). We have performed 
one simulation with 128 3 baryonic resolution elements, the same number of dark matter particles, 
and a number of stellar particles were formed during the simulation. The simulation box size was 
fixed to 3/i -1 Mpc, which resulted in the total mass resolution of 1.3 x 10 6 /i _1 Mq. For reference, 
the Jeans mass at z = 10 and T = 10 4 K is about 10 10 /i _1 Mq. The spatial resolution was fixed at 
1.5/i _1 comoving kiloparsecs. Because of the small box size, this simulation cannot be continued 
to z = 0. Instead, we stopped the simulation at z ~ 3. 

Even at this redshift our simulation is suffering from the lack of spatial of resolution, as 
can be illustrated by Figure [l], which shows the evolution of the central density for the dark 
matter, baryonic and stellar components of the most massive object from our simulation. The 
central density is defined as the average density within the sphere of two resolution lengths of 
our simulations (~ 900 pc at z = 3.8). As an illustration, we also show the central density of the 
Navarro, Frenk, and White (NFW) profile, which is defined in the same way and is therefore a 
function of resolution (and thus time, as we keep the resolution fixed in the comoving coordinates) 
- as redshift decreases, so does our spatial resolution in the physical coordinates. The fact that 
the NFW density actually starts to decrease for z < 6 indicates that our resolution becomes 
comparable to the characteristic radius of the object - the radius where the local density slope 
approaches -2. 

Another simulation with eight times more resolution elements (256 3 ) and 6/i _1 Mpc box was 
also performed. The large simulation thus had the same mass resolution as the small one, and 
the spatial resolution in the large simulation was fixed at \.2hr x comoving kiloparsecs. However, 
because the large simulation required a computational expense beyond what was available to us, 
it was terminated at z = 9.5. Thus, we used the large simulation to verify numerical convergence 
and estimate missing small scale power, but we will use the small (128 3 ) simulation as the source 
for scientific results. 

By comparing the large and small simulation, we have found that the small simulation 
included most of the small scale power that was initially present in the baryonic component. Thus, 
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Fig. 1. — Evolution of the central dark matter density (solid lines), baryonic density (dashed lines), and 
stellar density (dotted lines) for the most massive object in our simulation. The thin solid line shows the 
central total density of the same object if it had the Navarro-Frenk- White (NFW) density profile. The 
central core appears because of finite spatial resolution in simulations. 

our results are not significantly affected by the finite resolution in the initial conditions (/c-space 
resolution), but they are, of course, subject to finite mass and spatial resolution. 



3. Results 

Since we are concerned with the process of formation of galactic bulges and small ellipticals, we 
will focus in this paper on properties of individual objects formed in our simulations. Specifically, 
we will focus on four most massive objects. Each of those objects contains more than ten thousand 
particles of each kind (i.e. the dark matter, gas, and stars), and thus they are fully resolved 
numerically. 



- 5 - 

Table 1: Four Most Massive Objects at z = 4 



Object Total mass (M ) Baryonic mass ( M ) Stellar mass ( M ) M b /M t M*/M 6 



A 


4.0 x 10 10 


8.2 x 10 9 


1.5 x 10 9 


0.20 


0.18 


B 


2.9 x 10 10 


5.1 x 10 9 


1.2 x 10 9 


0.18 


0.23 


C 


1.8 x 10 10 


3.3 x 10 9 


0.9 x 10 9 


0.18 


0.26 


D 


1.4 x 10 10 


2.2 x 10 9 


0.6 x 10 9 


0.16 


0.29 



Table 2: Four Most Massive Objects at z = 3 



Object Total mass (M ) Baryonic mass ( Mq ) Stellar mass ( M ) M b /M t M*/M b 



A 


7.1 x 10 10 


1.1 x 10 10 


2.4 x 10 9 


0.15 


0.23 


B 


2.9 x 10 10 


4.5 x 10 9 


1.2 x 10 9 


0.16 


0.27 


C 


1.8 x 10 10 


3.0 x 10 9 


0.9 x 10 9 


0.17 


0.22 


D 


1.6 x 10 10 


2.5 x 10 9 


0.9 x 10 9 


0.16 


0.34 



Tables [T] and |2| present the general properties of the four objects: their total, baryonic, and 
stellar masses, as well as the ratio of the baryonic to the total mass, and the stellar to the baryonic 
mass at z ~ 4 and z ~ 3 respectively. 

A few observations can easily be made from the table. First, the objects contain more baryons 
than the cosmic average of Qb/^o = 0.13. Second, they are still dominated by gas, as only about 
20-30% of their baryonic mass is turned into stars when we terminated the computation at z ~ 3. 
These objects have a moderately high rate of star formation (several solar masses per year), and 
would appear to correspond well to the numerous small blue objects seen in the Hubble Deep 
Field (Contardo, Steinmetz, & Alvensleben 1998). 

But, it must be noted here that the amount of gas turned into stars really depends on the 
specifics of the star formation algorithm adopted in the simulation. We make no pretense that 
we can fully account for the star formation in our simulations, it is clear that modeling of star 
formation needs to be much advanced before simulations could claim to predict the star formation 
rate in (proto-)galaxies with sufficient detail. Thus, for the purpose of this paper we will pay 
little attention to the stellar component of our simulated bulges, and will focus primarily on the 
total baryonic component instead. We assume that with proper star formation algorithm, our 
simulated bulges will form stars at an appropriate rate as long as the total baryonic distribution 
is compatible with observations. 

Since at z < 3.8 our results become substantially contaminated by the lack of spatial 
resolution, we will restrict most of our analysis to the redshift z = 3.8, the lowest redshift at 
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Fig. 2. — Density profiles for the dark matter (solid lines), total baryons (gas and stars, dashed lines), and 
stars (dotted lines) of the four most massive objects at z = 3.8 as a function of radius in physical (not 
comoving) units. The bold vertical bar marks the spatial resolution of the simulation (450 pc). 

which characteristic radii of bound objects are still well resolved. Density profiles of the four most 
massive objects are shown in Figure |2[ The most massive object, object A, is less dense at the 
center than other three objects because it experienced a major merger shortly before z = 4 and has 
not fully relaxed yet. Object D has also experienced a major merger at z ~ 6, whereas objects B 
and C have been accreting matter quietly since z ~ 10. We point out here that in all four objects 
baryonic density at the center is greater than or similar to the dark matter density, i.e. baryons 
in all four objects are self-gravitating. As an aside we note that should the gaseous component 
dominate at the center at any time, a rapid collapse would occur which presumably would be 
accompanied by rapid star formation or, should a black hole of sufficient mass be presents, the 
concomitant flare up of an AGN. Clearly, we have not sufficient resolution to follow this phase 
(but see Abel, Bryan, & Norman 1998). 

Evolution of average properties of these objects is shown in Figure ||. One can immediately 
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Fig. 3. — Evolution of the total mass (solid lines), baryonic fraction (dashed lines), and stellar fraction (per 
unit baryonic mass, dotted lines) of the four most massive objects. 

see that all four objects are experiencing heavy merging at z ~ 4 — 6, increasing their mass by 
about an order of magnitude. The stellar fraction of object A decreased significantly at z = 4 — 5 
because of accretion of a large quantity of fresh gas. Also noticeable is the increase in the total 
baryonic fraction in object A at z ~ 4. Since baryons are almost self-gravitating at the center 
of object A, their efficiently cool and collapse toward the center, leading to the increase in the 
baryonic fraction of the object. 

Figure [l] shows the evolution of the central density for the dark matter, baryonic and stellar 
components of object A. One can see that the dark matter density does not change systematically 
with time (albeit fluctuating significantly) until the lack of resolution contaminates results, because 
the object is close to the virial equilibrium. On the contrary, the baryonic density increases with 
time because of efficient cooling at the center of the object. The recent merger at z ~ 5 triggered 
a considerable increase in the central density of gas, but this increase has not yet resulted in the 
burst of star formation. We expect, that if we continued the simulation to lower redshift with 



- 8 - 




R (kpc) R (kpc) 



Fig. 4. — Surface mass density profiles for the dark matter (solid lines) and total baryons (dashed lines) 
of the four most massive objects at z = 3.8 as a function of radius in physical (not comoving) units. The 
profiles are terminated at the resolution limit of the simulation (450 pc). Tilted long-dashed lines show r 1 / 4 
law for the baryonic profiles, and dotted lines show exponential profiles. Bold arrows show characteristic 
radii R e for r 1 / 4 profiles. The horizontal long-dashed line show the central surface density of a homogeneous 
top-hat sphere. 

higher spatial resolution, object A would experience a burst of star formation at the center, which 
would transform most of the gas into stars on a rather short time-scale. 

We are now ready to address the major question of this paper: are those objects formed in 
the simulation resemble real galactic bulges? In order to answer this question, we show in Fig. |I| 
the surface density profiles for our four objects. Also, for the total baryonic profile we compute 
the exponential fits in the form: 

E(R) = See 1 " R / Rc (3) 

and r 1//4 law fits, 

S(i?) = S e 10- 3 - 3307 (W fe ) 1/4 - 1 ). (4) 
As can be seen from Fig. ^| the central 3 kpc of the simulated bulges are equally well fitted both by 
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Table 3: Fit Parameters for the Four Most Massive Objects at z = 3.8 

Object R c (kpc) S c (M /pc 2 ) R e (kpc) S e (M /pc 2 ) 



A 0.70 1200 1.47 310 

B 0.54 1270 1.11 310 

C 0.58 664 1.27 140 

D 0.46 790 0.91 200 

Galactic bulge a 0.58 650 c 0.97 250 c 

Galactic bulge b 0.47 560 c 0.76 260° 



"Kent, Dame, & Fazio 1991 
6 Freudenreich 1998 

C A mass to light ratio of 3 is assumed in converting luminosity to mass. 

the exponential and by the de Vaucouleurs profiles, and we lack the resolution to distinguish them 
at the very center. We also note that our resolution is barely enough to resolve the characteristic 
scales of the two profiles, so we may overestimate the characteristic radii somewhat, but not by a 
large factor. However, both, Kent, Dame & Fazio (1991) and Freudenreich (1998) models for the 
galactic bulge are well fitted by the exponential profile for the range 0.1 kpc < r < 3 kpc (with 
respective rms errors of 3 and 8 percent respectively), and thus agree with our simulations over 
the range of scales which we can resolve. 

We put together the parameters of the fits in Table |3|. In addition, we list the parameters 
that describe the Galactic bulge from Kent, Dame, & Fazio (1991) and Freudenreich (1998) for 
comparison. As one can see, objects that we observe in the simulation are very similar to (or 
perhaps even a little bit larger than) the bulge of Milky Way, provided, they convert most of their 
gas into stars. 

Why did we then erroneously conclude in the Introduction that the CDM-type models predict 
too low density bulges? The answer to this puzzle is again illustrated in Fig. ^. The horizontal 
long-dashed line in that figure shows the central surface density for a homogeneous top-hat sphere, 
i.e. for a spherical object with the constant density of 200 times the average density of the universe 
and with a radius equal to the virial radius of object A (all four objects have quite similar virial 
radii). One can immediately see that the top-hat model underestimates the central density of an 
object by about two orders of magnitude! Even the density at the characteristic radii (Rc, Re) is 
30-100 times greater than the fiducial top-hat virial density. 

Evolution of the two fit parameters, the characteristic radius and density, is shown in Fig. 
H for the four most massive objects in the simulation. The densities increase steadily as the 
gas continues to accrete and cool inside the dark matter halos, whereas characteristic radii stay 
approximately constant for all objects. 
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Fig. 5. — Evolution of the characteristic radii (lower panel) and densities (upper panel) for the four objects: 
A (solid lines), B (long- dashed lines), C (short-dashed lines), and D (dotted lines). The horizontal shaded 
areas mark the range of values for the Galactic bulge, and the vertical line marks the boundary z = 3.8 
beyond which our simulation fails to resolve cores of the galactic bulges. 

Finally, if we want to demonstrate that we can form galactic bulges in a realistic cosmological 
simulation, we should address the question of the bulge shape. The bulge of our Galaxy is 
quasi-spherical, or, at the very least, slightly ellipsoidal. Is this shape also reproduced in the 
simulation? Figure || addresses this question. In it we show the axis ratios for the dark matter, 
gas, and stars for our four objects as a function of radius. One can see that in the central parts 
the gas and the stellar distributions are quite close to spherical. Shapes at larger radii, r > 3 kpc, 
shown shaded in Fig. ^| vary significantly among all objects, but at those large distances the gas 
is far from equilibrium, and the minute shape of its distribution has little relation to its final 
state. The fact that the objects we observe in the simulation are more-or-less spherical, rather 
than disk-shaped, is due to the fact that at high redshift the slope of the linear power spectrum 
of the density fluctuations is close to —3. This results in a range of scales becoming nonlinear 
almost at the same time, which, in turn, leads to heavy merging among objects observed in the 



-11- 




Fig. 6. — Axis ratios for the dark matter (solid lines), total baryons (gas and stars, dashed lines), and stars 
(dotted lines) of the four most massive objects at z = 3.8 as a function of radius in physical (not comoving) 
units. 

simulation. Thus, gaseous disks do not have enough time to form, and the shape of the objects 
remains quasi-spherical. Only at later times, when the rate of merging falls down, can a gaseous 
disk form inside an object. 

In order to confirm this claim, we show in Figure ^ the velocity dispersion and temperature 
profiles for the four objects mentioned above. As one can see, the gas velocity dispersion is 
somewhat larger than the gas temperature, and so the turbulent motions in the gas provide more 
than 50% of the support against the gravity, because the molecular and/or magnetic viscosity is 
not large on the scales which we can resolve, it would take some time before the gas can settle 
into a rotationally supported disk. Thus, in order for our objects to resemble stellar bulges at the 
present epoch, the turbulence dumping time should be large compared to the star formation time. 
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Fig. 7. — Velocity dispersion and gas temperature profiles for the four most massive objects at z = 3.8. 
Shown are the dark matter velocity dispersion (solid lines), gas temperature (long-dashed lines)v, gas velocity 
dispersion (short-dashed lines), and stellar velocity dispersion (dotted lines). The velocity dispersions are 
converted into temperature units. 

4. Conclusions 

We have showed that objects that form in a realistic cosmological simulation of a CDM-type 
cosmological model do look similar to bulges of normal galaxies. Our objects are still 75% gaseous 
at z ~ 4, but they form stars at a high rate, and when most of the gas in those objects will be 
converted into stars by redshift z = 2. In this time frame they resemble the numerous small blue 
galaxies seen in the HDF. The number density of such objects predicted from our simulation is 
about a factor of two larger than the actual number observed in the HDF and LBG galaxy samples 
(Steidel et al. 1999), however the precise comparison with the observations depends crucially 
on the (highly uncertain) assumptions about the star formation, and is not possible at the this 
moment. Later the simulated objects will look like slightly ellipsoidal stellar objects with the 
density profiles well fit by the exponential profile and with the parameters of the fit similar to the 
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parameters of the Milky Way bulge and small elliptical galaxies. 

We thus conclude that currently favorable CDM-type cosmological models have no difficulty 
in reproducing observed properties of galactic bulges. On the contrary, models that have galaxy 
formation at z ~ 30 (Peebles 1997) would form bulges that are two to three orders of magnitude 
more dense than the observed ones. This conclusion is further boosted by the consideration that, 
due to the limited mass and spatial resolution of our simulation, we can only underpredict the 
densities and masses of cosmological objects. With several bulges observed in approximately 
correct mass range (note the local group with M31, M33, and the Galaxy), and a total luminosity 
of about L* in the volume of 27/i~ 3 Mpc 3 , we conclude that the model produces approximately 
the correct density of bulges. 

We are grateful to the referee Andreas Burkert for the valuable comments that significantly 
improved the paper. This work was supported in part by the UC Berkeley grant 1-443839-07427 
and the NCSA/NSF Subaward 766/ASC97-40300. Simulations were performed on the NCSA 
Origin2000 supercomputer under the grant AST-970006N. 
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